// @HEADER
// ***********************************************************************
///
//             SiYuan: A numerical PDE solver
//                 Copyright (2022) YUAN Xi
//
// ***********************************************************************
// @HEADER

#include "Teuchos_ParameterList.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_Assert.hpp"
#include "Phalanx_DataLayout_MDALayout.hpp"
#include "Phalanx_FieldManager.hpp"
#include "Panzer_PhysicsBlock.hpp"

#include "Panzer_BasisIRLayout.hpp"

// Evaluators
#include "Panzer_ConstantFlux.hpp"
#include "Panzer_Integrator_BasisTimesScalar.hpp"
#include "Panzer_Integrator_BasisTimesVector.hpp"
#include "Panzer_Constant.hpp"
#include "Convection.hpp"
#include "Impedance.hpp"

#include "Phalanx_MDField.hpp"
#include "Phalanx_DataLayout.hpp"
#include "Phalanx_DataLayout_MDALayout.hpp"

#include <bits/stdc++.h>

// ***********************************************************************
template <typename EvalT>
SiYuan::BCStrategy_Neumann<EvalT>::
BCStrategy_Neumann(const panzer::BC& bc, const Teuchos::RCP<panzer::GlobalData>& global_data) :
  panzer::BCStrategy_Neumann_DefaultImpl<EvalT>(bc,global_data), m_gd(global_data)
{}

// ***********************************************************************
template <typename EvalT>
void SiYuan::BCStrategy_Neumann<EvalT>::
setup(const panzer::PhysicsBlock& side_pb,
      const Teuchos::ParameterList& /* user_data */)
{
  const std::string dof_name = this->m_bc.equationSetName();
  // need the dof value to form the residual
  //this->requireDOFGather(dof_name);

  const std::string residual_name = "Residual_" + this->m_bc.identifier();
  const std::string flux_name = "Neumann_" + this->m_bc.equationSetName();
  const std::map<int,Teuchos::RCP< panzer::IntegrationRule > >& ir = side_pb.getIntegrationRules();
  TEUCHOS_ASSERT(ir.size() == 1); 
  const int integration_order = ir.begin()->second->order();

  if( this->m_bc.strategy() == "Acoustics Velocity" || this->m_bc.strategy() == "Impedance" ) 
  {
    std::string resr_name = residual_name + "_Real";
    std::string dofr_name = dof_name + "_Real";
    std::string fluxr_name = flux_name + "_Real";
    this->addResidualContribution(resr_name,dofr_name,fluxr_name,integration_order,side_pb);

    std::string resi_name = residual_name + "_Img";
    std::string dofi_name = dof_name + "_Img";
    std::string fluxi_name = flux_name + "_Img";
    this->addResidualContribution(resi_name,dofi_name,fluxi_name,integration_order,side_pb);
  }
  else
    this->addResidualContribution(residual_name,dof_name,flux_name,integration_order,side_pb);
}

// ***********************************************************************
template <typename EvalT>
void SiYuan::BCStrategy_Neumann<EvalT>::
buildAndRegisterEvaluators(PHX::FieldManager<panzer::Traits>& fm,
			   const panzer::PhysicsBlock& pb,
			   const panzer::ClosureModelFactory_TemplateManager<panzer::Traits>& /* factory */,
			   const Teuchos::ParameterList& /* models */,
			   const Teuchos::ParameterList& /* user_data */) const
{
  const std::vector<std::tuple<std::string,std::string,std::string,int,Teuchos::RCP<panzer::PureBasis>,Teuchos::RCP<panzer::IntegrationRule> > > 
    data = this->getResidualContributionData();

  std::string residual_name = std::get<0>(data[0]);
  std::string dof_name = std::get<1>(data[0]);
  std::string flux_name = std::get<2>(data[0]);
  Teuchos::RCP<panzer::IntegrationRule> ir = std::get<5>(data[0]);
  Teuchos::RCP<const panzer::FieldLayoutLibrary> fll = pb.getFieldLibrary()->buildFieldLayoutLibrary(*ir);
  Teuchos::RCP<panzer::BasisIRLayout> basis = fll->lookupLayout(dof_name);
  Teuchos::RCP<const panzer::PureBasis> pbasis = fll->lookupBasis(dof_name);
//  bool aa = pbasis->isVectorBasis();

  // provide a constant flux target value to map into residual
  Teuchos::ParameterList p("Neumann BC");
  bool is_complex = false;

  Teuchos::RCP< PHX::Evaluator<panzer::Traits> > op;
  try {
    if( this->m_bc.strategy() == "Flux" ) {
      p.set("Name", flux_name);
      if( basis->getBasis()->isVectorBasis() )
      {
        p.set("Data Layout", ir->dl_vector);
        p.set("Value X", this->m_bc.params()->template get<double>("Value X"));
        if(pbasis->dimension()>1)
          p.set("Value Y", this->m_bc.params()->template get<double>("Value Y"));
        if(pbasis->dimension()>2)
          p.set("Value Z", this->m_bc.params()->template get<double>("Value Z"));
    
        op = Teuchos::rcp(new panzer::ConstantVector<EvalT,panzer::Traits>(p));
      }
      else {
        p.set("Data Layout", ir->dl_scalar);
        p.set("Value",  this->m_bc.params()->template get<double>("Value"));
        op = Teuchos::rcp(new panzer::Constant<EvalT,panzer::Traits>(p));
      }
    }
    else if( this->m_bc.strategy() == "Convection" )
    {
      p.set("Name", flux_name);
      p.set("Ambient Temperature", this->m_bc.params()->template get<double>("Ambient Temperature"));
      p.set("Heat Transfer Coefficient", this->m_bc.params()->template get<double>("Heat Transfer Coefficient"));
      p.set("DOF Name", dof_name);
      p.set("IR", ir);
      p.set("Basis", basis); 
      op = Teuchos::rcp(new Convection<EvalT,panzer::Traits>(p));
    }
    else if( this->m_bc.strategy() == "Radiate" )
    {
      p.set("Name", flux_name);
      p.set("Ambient Temperature", this->m_bc.params()->template get<double>("Ambient Temperature"));
      p.set("Radiate Coefficient", this->m_bc.params()->template get<double>("Radiate Coefficient"));
      p.set("DOF Name", dof_name);
      p.set("IR", ir);
      p.set("Basis", basis); 
      op = Teuchos::rcp(new Radiate<EvalT,panzer::Traits>(p));
    }
    else if( this->m_bc.strategy() == "Impedance" )
    {
     // auto gd = this->getGlobalData();
      double omega;
      try {
        omega = m_gd->constants.at("Frequency");
      }
      catch(std::out_of_range&) {
        std::cout << "Definition of Frequency MUST be given!" << std::endl;
      }
      is_complex = true;
      flux_name.erase(flux_name.end()-5,flux_name.end());
      p.set("Name", flux_name);
      auto vn  = this->m_bc.params()->template get<Teuchos::Array<double>>("Impedance");
      p.set("Impedance", vn);
      p.set("Density", this->m_bc.params()->template get<double>("Density"));
      p.set("Sound Speed", this->m_bc.params()->template get<double>("Speed of Sound"));
      p.set("DOF Name", dof_name);
      p.set("IR", ir);
      p.set("Basis", basis); 
      omega = 2.0*M_PI*omega;
      p.set("Angular Frequency", omega);
      op = Teuchos::rcp(new AcousticImpedance<EvalT,panzer::Traits>(p));
    }
    else if( this->m_bc.strategy() == "Acoustics Velocity" )
    {
    //  auto gd = this->getGlobalData();
      double omega;
      try {
        omega = m_gd->constants.at("Frequency");
      }
      catch(std::out_of_range&) {
        std::cout << "Definition of Frequency MUST be given!" << std::endl;
      }
      is_complex = true;
      flux_name.erase(flux_name.end()-5,flux_name.end());
      p.set("Name", flux_name);
      p.set("Density", this->m_bc.params()->template get<double>("Density"));
      p.set("Sound Speed", this->m_bc.params()->template get<double>("Speed of Sound"));
      auto vn  = this->m_bc.params()->template get<Teuchos::Array<double>>("Normal Speed");
      p.set("Normal Speed", vn);
      p.set("DOF Name", dof_name);
      p.set("IR", ir);
      p.set("Basis", basis); 
      omega = 2.0*M_PI*omega;
      p.set("Angular Frequency", omega);
      op = Teuchos::rcp(new AcousticVelocity<EvalT,panzer::Traits>(p));
    }
    else
      throw std::runtime_error("Neumann Strategy Not Defined!");
  }
  catch (std::exception& e) {
      std::cout << e.what() << std::endl;
  }
  this->template registerEvaluator<EvalT>(fm, op);

  // add contribution to the residual 
  double multiplier(1);
  if( is_complex )    // we need set a complex variable here
  {
    residual_name = std::get<0>(data[0]);
    flux_name = std::get<2>(data[0]);
    Teuchos::RCP<PHX::Evaluator<panzer::Traits>> opr = Teuchos::rcp(new
        panzer::Integrator_BasisTimesScalar<EvalT, panzer::Traits>(panzer::EvaluatorStyle::EVALUATES,
        residual_name, flux_name, *basis, *ir, multiplier));
    this->template registerEvaluator<EvalT>(fm, opr);

    residual_name = std::get<0>(data[1]);
    flux_name = std::get<2>(data[1]);
    Teuchos::RCP<PHX::Evaluator<panzer::Traits>> opi = Teuchos::rcp(new
        panzer::Integrator_BasisTimesScalar<EvalT, panzer::Traits>(panzer::EvaluatorStyle::EVALUATES,
        residual_name, flux_name, *basis, *ir, multiplier));
    this->template registerEvaluator<EvalT>(fm, opi);
  } 
  else
  {
    if( basis->getBasis()->isVectorBasis() ) {
      Teuchos::RCP<PHX::Evaluator<panzer::Traits>> op = Teuchos::rcp(new
      panzer::Integrator_BasisTimesVector<EvalT, panzer::Traits>(panzer::EvaluatorStyle::EVALUATES,
      residual_name, flux_name, *basis, *ir, multiplier));
      this->template registerEvaluator<EvalT>(fm, op);
    }
    else {
      Teuchos::RCP<PHX::Evaluator<panzer::Traits>> op = Teuchos::rcp(new
        panzer::Integrator_BasisTimesScalar<EvalT, panzer::Traits>(panzer::EvaluatorStyle::EVALUATES,
        residual_name, flux_name, *basis, *ir, multiplier));
      this->template registerEvaluator<EvalT>(fm, op);
    }
    
  }

}

// ***********************************************************************
template <typename EvalT>
void SiYuan::BCStrategy_Neumann<EvalT>::
postRegistrationSetup(typename panzer::Traits::SetupData /* d */,
		      PHX::FieldManager<panzer::Traits>& /* vm */)
{
  
}


// ***********************************************************************
template <typename EvalT>
void SiYuan::BCStrategy_Neumann<EvalT>::
evaluateFields(typename panzer::Traits::EvalData /* d */)
{
  
}
